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In this paper we study synchronized motions in complex networks in which there are distinct 



ferent for nodes in different groups. Both continuous time and discrete time systems are considered. 
We initially focus on the case where two groups are present and the network has bipartite topology 
' (i-e., links exist between nodes in different groups but not between nodes in the same group). We 

also show that group synchronous motions are compatible with more general network topologies, 
, — where there are also connections within the groups. 



I. INTRODUCTION 



Because of its relevance in a wide variety of physical, biological, social and engineering contexts, synchronization of 



^ ' complex networks of coupled dynamical systems has recently received increasing attention. In this paper, we analyze 
synchronized (possibly chaotic) motions in complex networks of coupled groups of dynamical systems. Here by a 
group we mean a collection of systems that have the same dynamics, with any given group consisting of systems with 
dynamics that is different from the dynamics of systems in the other groups. Specifically, we will show that under 
certain circumstances, multiple group-synchronous evolutions may exist in such networks. In this type of synchronous 
motion, the evolution of the states of systems within a particular group are the same, while the states of members of 



d ' different groups although coherently related, are in general different (indeed the state vectors of systems in different 
groups may have different dimensionality). 

The problem of collective behavior in a network connecting members of different groups is of broad interest. As a first 
example, we note that many efforts have been devoted to the study of teams (groups) of interacting robots performing 
synchronous coordinated tasks 1 , . Other studies have regarded the coordination and control of several squadrons 
(groups) of unmanned autonomous vehicles to accomplish interdependent tasks, such as cooperative searches and 
attacks 0,0,01. In the social networks literature, distinct collective behaviors of individuals are often related to their 



sex, social status and/or race. Some studies have clearly pointed out how men's and women's social behaviors differ, 
even in situations where they are found to interact tightly, as in virtual communities or internet chats. In the brain 
functional assemblies of neurons have been observed to display distinct interdependent synchronous oscillations J] . 
Collective dynamics of groups displaying multi-synchronous behaviors have also been uncovered in ecological systems 
7|, lalSl, where competition could have favored the evolution of different synchronous behaviors of different species. 



For example, some corals are known to spawn synchronously during a particular season of the year 



10|. At the same 



time, different coral species typically spawn in different months, possibly to prevent hybridization of the species and/or 



Distinct roles of males and females (as in the 



as a mechanism to relieve larvae from interspecific competition 
case of social networks) infiuence the sexual activity of animals, where reproductive synchrony has been speculated 
to benefit survival of progeny by decreasing the likelihood of the male deserting his partner. 

Phase synchronization between essentially different chaotic systems has been the subject of intensive study since 
the appearance of the paper 



12l |. However, here we will be interested in complete (full) synchronization. Moreover, 
multiple synchronized motions of identical oscillators have been observed to coexist in complex networks characterized 
by strong community structure 
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15|. Here we will show that under certain conditions, multiple synchronous 



behaviors of systems with group properties can occur 



24|. 



In this paper we focus on the case where there are two groups. In Sec. II we consider bipartite network topology 
and continuous time dynamics and present examples of both periodic and chaotic synchronous behavior. In Sec. HI, 
discrete time systems are discussed. In the case of bipartite network topology studied in Sees. II and III, a compact 
master stability function 



16[ description for evaluating the stability of group synchronous motions is possible. When 
there are more than two groups or when there are two groups but a non-bipartite network structure, the stability 
analysis is generally more difficult. In Sec. IV, we remove the constraint of bipartite network topology, and we show 
that the stability of the multi-synchronous evolutions is indeed possible under these more general conditions and that 
it can be enhanced when connections are allowed between systems belonging to the same group. 



II. CONTINUOUS TIME BIPARTITE SYSTEMS 

In this section we focus on continuous time systems and consider a bipartite network connecting two groups. We find 
the conditions that allow a synchronization manifold and study its stability by means of a master stability function 
approach. 
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A. 



Formulation 



The individual equation of an isolated (uncoupled) node is denoted by ii = F(xi), i = 1, N^, for the nodes in the 
first group Sx and by yj = G{yj) j — 1, Ny for the nodes in the second group Sy, where Xi (yj) is an n3;-diniensional 
(nj,-dimensional) state vector and F : H"" — > R"^ and G : R"" ^ R"". The dynamical equations of the network 
systems are as follows: 



where A is an iV^; x Ny coupling matrix, whose entries {^ij} represent the intensity of the direct interaction from 
system j in Sy to i in Sx- Analogously the entries {Bji} of the Ny x Nx matrix B represent the interaction from 
system i in Sx to j in Sy. The interaction function H [L) is a mapping from R"" to R""= (from R"=" to R"" ). 

We now consider the possibility of the existence of multi-synchronous solutions, where by this we mean that 
xi{t) — X2{t) ~ ... = XN^{t) = Xs(t) and yi(t) ~ j/2(i) — ■■■ — UNyit) = ys{t)- Substituting such an assumed solution 
in ([1]) , we see that in order for a multi-synchronous state to exist the sum Aij must be independent of i and 
the sum Bji must be independent of j. If we denote the first sum by a and the second sum by b, then by the 
replacements aH H and A/a^ A {bL L and B/b B) we see that, without loss of generality, it suffices to set 
a = 6 = 1, 




(1) 



% = G{yj) + Bj,L{x,), J = 1, iV, 



i=l 



AT, 



y 




(2a) 




(2b) 



i=l 



Thus the equations of motion for the synchronized dynamics are 



Xs^F{xs)+H{ys), 



(3) 



y's = G[ys) + L{xs)- 
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B. Synchronization Stability 

In what follows we seek to characterize the stability of the above defined synchronous state. Linearization of the 
system ([Ij around the synchronous evolutions Xs{t) and yields: 



dx, = DF{xs)5x, + ^A,jDH{y,)5yj, i = l,...,N^ 



(4) 



Syj = DG{y,)5yj + ^ B.j,DL{xs)5x,, j = 1, ...,Ny. 



The Lyapunov exponents of the dynamics of a synchronous state (xs{t), ys{t)) are those associated with the follow 
system: 

6xs = DF{xs)Sxs + DH{ys)6ys, 

(5) 

5ys = DG{ys)5ys + DL{xs)Sxs, 

obtained by linearization of Eqs. Note that the synchronous evolutions Xs and ys might, e.g., be stationary, 

periodic, or chaotic. 

We now assume that the (N^ + Ny) independent solutions of (jH) can be expressed in the form Sxi — Cri-.Sx, 
i = l,...,Nx and Syj = Cy^Sy, j — l,...,Ny, where {c^^.} and {cy^} are appropriate time-independent scalars. This 
assumed form will encompass all possible linear solutions of (4) if the space of vectors given by the possible values of 
Cxi,Cy^{i — 1, ...,Nx;j — 1, ...,Ny) has dimension + Ny. As we shall see, this is the case (cf. Eq. (9) to follow). 
With the assumption, 6xi = CxiSx, i = 1, ■■■,Nx and Syj = Cy-6y, j = 1, ■■■,Ny, Eqs. ^ become 

Ny 



Cx.Sx = Cx,DF{xs)6x + (^AijCy^)DH{ys)Sy, i^l,...,Nx (6a) 
Cy^Sy = Cy^DG{ys)6y + (^B.jiCx^)DL{xs)6x, j^l,...,Ny. (6b) 



Thus in order that ([5a|) . (respectively ([6b)) ). is satisfied for all i, (respectively j), we require that Cx^J2j ^ij^Vj — 
where v is independent of i and Cy^'^^ BjiCx^ = i], where 77 is independent of j. After defining the vectors = 
(cxi,Cx2, ■■■^Cxj^J'^ and 

% ~ (cyi J : Cj^jv )^: these conditions may be rewritten as: Acy — vcx and Bcx — T/Cy; 

that is, 
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Using this in (5) we obtain 




(7) 



Sx = DF{xs)6x + vDH{ys)5y, 
Sy ^ DG{ys)Sy + t]DL{xs)Sx. 



(8) 



One particular solution of ^ is obtained when i/ = rj = \, i.e., 



Q 




(9) 



where A belongs to the set A = {A^}, i — I, ---^Nx + Ny of the (possibly complex) eigenvalues of the matrix Q. 
Rewriting © as 





(10) 



shows that solution of ^ yields all the possible solutions of ((7]) by setting v = \z , -q ^ \/ z ^ — c^^, Cy — zCy, where z 
is a free parameter. Also since A and B are real, the spectrum of Q is symmetric about the i?e(A) axis. Furthermore, 
we note that, if A is an eigenvalue of Q, then, by letting z — —1, we see that —A is also an eigenvalue. Thus the 
spectrum of Q is symmetric about both the Re{X) axis as well as the Im{\) axis. 

Moreover, the stability of the synchronous evolutions docs not depend on the particular z. In fact, if in Eqs. ^ 
we let V = \z, rj = A/z, 6y — zSy, we see that Sx, Sy satisfy Eqs. ([8]) with v — rj = X. Thus it suffices to consider ([9|), 
and we rewrite Eqs. ^ as 
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5x = DFix,)5x + XDH{ys)Sy, 

(11) 

S^^DG{y,)Sy + XDL{xs)Sx, 

n 

where A — Ai, A2, Ajy. We define a master stability function 16] for this problem, denoted M{Xi), where M 
associates to A, the maximum Lyapunov exponent of the system pip . Note that the function A/(A) can be determined 
without knowledge of the matrix Q. Thus the synchronization stability problem is decomposed in two parts, (i) a part 
dependent only on the couplings H and L and on the individual system dynamics F and G, but not on the network 
topology (i.e. not on the matrix Q), and (ii) a part dependent solely on the network topology (determination of the 
spectrum of Q). 

Another important consequence of the invariance of ([8|) under the transformation 1/ ^ Xz, rj ^ X/ z is that the 
synchronous state stability for an eigenvalue A is the same as for —A {z ^z). Thus only those eigenvalues with, 
e.g., Re{X) > need to be tested. 

C. Spectrum of Q 

The matrix Q has a pair of real eigenvalues 1 and — 1 . This follows because the sums of the components for all rows 
of A and B are one. The eigenvalue +1 corresponds to an eigenvector all of whose components have the same value; 
while the eigenvector —1 corresponds to an eigenvector whose first components have the same value and whose 
remaining Ny components all have the negative of this value. Hence the eigenvalues +1 and —1 are associated with the 
directions parallel to the synchronization manifold; thus they may result in positive Lyapunov exponents corresponding 
to chaotic dynamics taking place in the synchronization manifold (xi — X2 — ... — XN^,yi =2/2 = ••• = VNy)- In order 
to check the stability of the synchronous evolutions, one should evaluate the master stability function M{Xi) for the 
remaining + Ny — 2 eigenvalues of Q, representing the stability of the motions transverse to the synchronization 
manifold. The synchronized state is stable if, M{Xi) < for all Xi in the set A' = {A — { — 1, +1}}- 

From the fact that the sum of the elements in every row of Q is one, with all zero elements on the main diagonal, the 
Gershgorin circle theorem implies that the spectrum of Q lies in the disc of unit radius in the complex plane, having 
its center at (0, 0). Note also that the matrix Q has at least \Nx — Ny\ zero eigenvalues, so that zero eigenvalues must 
always occur unless iV^; = Ny. In particular, if iV^ ^ Ny, a necessary condition for the stability of the synchronized 
coupled systems is that the Lyapunov exponents resulting from dx = DF{xs)Sx and dy = DG{ys)5y are all negative. 
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FIG. 1: A simple 5-nodes network. 

(Note that these exponents depend on H and L because the synchronous time evolutions Xs{t) and ys{t) depend on 
H and L.) In order to see that Q has at least \Nx — Ny\ zero eigenvalues, assume that > Ny. Then the rows 
of A (each of which has Ny < components) can span a space of at most dimension Ny. Hence the N^ + Ny = N 
rows of Q can span a space of most dimension 2Ny, and there are at least {N^ + Ny) — 2Ny = N^ — Ny independent 
homogeneous linear relationships between the rows of Q, implying that there are at least N^ — Ny zero eigenvalues. 

Moreover, the spectrum of Q can be obtained through the computation of the eigenvalues of the lower dimensional 
of the two matrices, AB and BA. In fact, by noticing that is a block diagonal matrix of the form 




(12) 



we have that, if A is in the spectrum of Q, then must be one of the A^^ eigenvalues of AB and/or one of the Ny 

eigenvalues of BA. Say A^min = niin i^x, Ny), define the Nmin x Nmin matrix, 



D 



AB, if Nx<Ny, 
BA, if Ny<N^, 



(13) 



and denote the spectrum of 2? by A = {Ai, Ajv^^^}. Then, since the eigenvalues of Q"^ are the square of the 
eigenvalues of Q, we have that the spectrum of Q is 



A = [0, 0, 0] U[± v^, ±^2, ± VV~] 



(14) 



where [0, 0, 0] denotes | A^^ — Ny\ zeros. (Note that by Eqs. (2) one of the eigenvalues of D is +1, corresponding to 
an eigenvector [1, 1, 1]"^.) 
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As an example, we consider the network in Fig. [T] For this network = 2, Ny = 3, and 



A 




B 



1 I 

2 2 

yo I J 



(15) 



Note that the row sums of A and B are one as required by Eqs. (|2al2bp . Since — 2 < Ny = 3, = AB^ so that 



D 




(16) 



The eigenvalues of this 2x2 matrix are 1 and w. Thus since \Nx ~ Ny \ — 1. Eq. yields the real spectrum, 



A = [-1, -\/w,0, ^/uJ, 1] 



(17) 



We now consider the spectrum of Q for large networks of two types: (i) Q is random, and (ii) Q is constrained to 
have a real spectrum but is otherwise random. To construct the matrix Q in these two cases we start with a matrix 
Q' of the form 



Q' = 



A' 



B' 



(18) 



and then take Q to be 



(19) 



where is a diagonal matrix with Ka ~ Q'^ . Equation (|19p insures that the row sums of Q are one as required 
by Eqs. (|2al2bp . For case (i) we choose the elements of A {B') randomly to be one with probabilities p^y (Pyx), 
and zero otherwise. Note that, in case (i) there is no correlation between Aij and Bji. For case (ii) we choose A' 
randomly with A^j = 1 with probability p and A'^j = otherwise, and we then set B' = {A')'^ . Thus in this case, Q' is 
symmetric and the first elements Ku of K arc the row sums of A' , while the next Ny components are the column 
sums of A' . For case (ii) the spectrum is real, since by multiplying by K^^^ , the eigenvalue equation Qc = Ac can be 
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rewritten as (K'^I'^Q' K'^l'^^d = Ac', where c' = K^^'^c. Because K~^I'^Q' K'^l"^ is symmetric, if Q' is, we see that 
the eigenvalues of Q are real in case (ii). Next we use numerical experiments to investigate the general properties of 
the spectrum of large random matrices Q in the above two cases. 

First we consider case (i). We take — Ny — 500 and find the spectrum of Q for randomly generated matrices 
with several values oi Pxy and pyx- Results are shown in Fig. [2l We see that there are two eigenvalues at A = ±1 and 
that all the other eigenvalues lie within a circle whose radius decreases as the average node degree increases (i.e., as 
Pxy and Pyx increase). We evaluated the scaling of Xmax — max^ |Ai| for Xi e A' with the network size N in the simple 
case where Nx — Ny = N/2 and pxy = Pyx = P- We hypothesize a scaling of the form Xmax{N,pxy) — CN^ , and 
perform numerical simulations with p ranging between and 1, and N ranging between 200 and 2000. Our numerics 
show that e ~ 1/2 (we note, however, that for pxy — > 1, Xmax = 0, independent of N). Thus, by assuming a scaling 
of the form Xmax{N ^p) = CN^^I"^ we obtained different values for C, as function of the probability p. In Fig. [3] the 
values of the logarithm of Xmax I C are shown to collapse to a straight line of slope —1/2 for different values of p, when 
plotted versus the logarithm of the network dimension N . The inset of FigO shows C versus p for p ranging between 
and 0.9 in steps of 0.1. The scaling Xmax ^ N^^^^ implies that, with increasing N, the spectrum A' shrinks toward 
the point (0,0). Moreover, A' also shrinks toward (0,0) as Pxy and Pyx approach one, independently of the network 
dimension A^. Thus both in the case of a very large network (i.e., A^ large) or complete network (i.e., Pxy^Pyx ^ 1), 
the whole spectrum of the eigenvalues in A' collapses onto the real eigenvalue 0. 

We now consider case (ii) where the spectrum of Q is real. Analogous to the results in Figs. [2]and[31 we find that for 
large A'^ the eigenvalues of Q in A' are distributed along the real line lying in a symmetric range, —Xmax l£ X < Xmax, 
where Xmax decreases toward zero with increasing p (see Fig. SJa)) as well as increasing A^ (see Fig. Hl^b)), with 
Xmax ~ N-^/^ for large N. 

The result that Xmax decreases with A^ and p (or pxy and Pyx) for both cases (i) and (ii) is quite significant. 
In particular, if Xmax ^ 1, then, for the purposes of evaluating the master stability function, it becomes a good 
approximation to set A = 0. This is a great simplification in that the master stability function now need be evaluated 
only for a single value of A, and its determination reduces to a computation on two uncoupled systems, 

Sx — DF(xs)Sx, 

(20) 

= DG{ys)Sy, 

where we again emphasize that, although the coupling functions H and L do not appear explicitly in (|20p . M still 
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Re{^.} 

FIG. 2: (Color online) Random networks with A'^; = Ny — 5 x 10^ nodes. The spectrum of Q is shown for three sets of values 
of Pxy and Pyx- Yellow (light gray) is used for p^y — Pyx = 0.05, red (dark gray) for pxy ~ 0.5 and pyx = 0.05, black for 
Pxy ~ Pyx ~ 0.5. The continuous line is used to represents the circle of unit radius centered at (0,0). The eigenvalues A — ±1 
associated with perturbations in the synchronization manifold are denoted by solid black dots. 

depends on H and L because the synchronous time evolutions, Xs{t) and ys{t), depend on H and L (see Eq. ([3]) ). 

D. Examples 

Example 1: Synchronized Periodic Motion. 

We consider the following coupled network dynamical equations, that are in the form ([1]), 

Ny 

ij(2) = -a;j(i) -a;j(2)(a;-(i) +a;-(2) - 1), i = l,...,7V^. (21) 
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FIG. 3: (Color online) log-log plot of Xmax/C versus N, for different values of p^y ranging between 0.1 and 0.9 in steps of 0.1. 
The straight line has slope —1/2. The inset shows C versus p. 



Vjil) = yj(2) + CTy X! 

1=1 

y,(2) = - 0.22/,(2)(2/|(i) - 1), j = 1, -.Ny. (22) 



In the absence of coupling (Tx = cry = 0, Eqs. (|2ip and (|22p both individually have global attractors on which the 
motion is periodic (i.e., they are limit cycle attractors). In particular, with ay — 0, Eq. (|22p is the Van der Pol 
equation. 

In order to measure the extent to which synchronization is achieved, we have monitored the asymptotic time average 
of the following two quantities: = i% Eil'i Yfj=ii\^i(i) + N»(2) -Xj(2)\) and Ey = -^Yh-i Ejl^i(ly»(i) " 

Vjil) I + \yi{2) ^ yj(2) 1); as functions of the control parameter Ux with Oy — 0.65. For each i and j we consider randomly 
chosen initial conditions in |a;ii_2| < 3 and |yji,2| < 3 and evolve the system for a long time (from t = Q to t — 300). 

The case of a network with a real spectrum, obtained as explained in Sec. II. C for the case (ii), is shown in Fig. 
[51 For this network we take = 200, Ny = 300 and p = 0.5, for which we find Xmax = 0.13. The upper panel of 
Fig. O shows Ex + Ey as functions of Gx for ay — 0.65 at different simulation times t = 100, 200, 300. We see that for 
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FIG. 4: The behavior of Xmax as function of p and of A'' for large random networks with real spectra (Q = K ^Q' and Q' 
symmetric), (a) A max versus for Nx ^ Ny ^ 5 X 10^ and Ny = 5 x 10^ nodes, (b) Log- log plot of Xmax/C versus A'^ for 
p = 0.1 and p = 0.5 showing that (as for case (i)) Xmax scales as N~^^'^ for large A*' (the solid line has slope —1/2). 

< fJa; < 0.4 the error decreases with time to very low values, indicating stable synchronization in this range of 
(the synchronized motion in this range is observed to be periodic). The lower panel shows the corresponding master 
stability function evaluated at A = (continuous line), and at A = Xmax = 0.13 (dashed line). We observe that the 
value for the zero crossing of M(A) is approximately at 0.4, for both A = and A = Xmax- For other values of A in the 
range < A < Xmax the curves are similar and have ax values at the zero crossings of the master stability function at 
approximately 0.4. Thus we find that the master stability function (lower panel of Fig. [5]) predicts a stable ax range 
of < da; < 0.4 in excellent agreement with our results from the full nonlinear computation (upper panel of Fig. [5|). 
Example 2: Synchronized Chaotic Motion. 

We now consider the following coupled network dynamical equations, 



Xi = -r{xi + h{xi)) 



i = l,...,Nx, 



(23) 



i=l 

yj{2) = -ivjii), 



(24) 



where h{x) = mix 



-{\x + 1| - Ix - 1|) and we take r = 4.6, q = 6.02, mo = -8/7, toi = -5/7. 



When Xi = Xs^i, yj(i.2) — ys(i,2)Vj, the three equation system formed by ([23|) and ([24|) has three attractors; two 
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FIG. 5: The upper panel shows + Ey versus a:c for (jy = 0.65. The continuous thin hne represents Ex + Ey a,t t = 100; 
the dashed hne, E^ + Ey a,t t = 200; the thick continuous hne, Ex + Ey at t = 300. The network parameters are as follows: 
Nx = 200 and Ny = 300, and p — 0.5. The lower panel shows the master stability fonction evaluated at A = (continuous 
line), and at A = Xmax = 0.13 (dashed line) versus ax- 

n 

are stable fixed points at (a;, y(i), 2/(2)) = (±3/2, 0, =f3/2) and the third is a chaotic attractor 17]. Thus, depending 
on the initial conditions, the motion in the synchronization manifold can be chaotic. We now investigate the stability 
of the synchronous chaotic motions for large N. For N >> 1 all the eigenvalues in A' tend to 0, and the synchronous 
evolution thus is stable if the maximum Lyapunov exponents associated with the following two (uncoupled) systems, 



Sjs ~ K{t)6x, where K{t) = — r — 



rmi, if |a;s| > 1, 



(25) 



and 



^y{i) = -'52/(1) +Sy{2), 



(26) 



are both negative. Note that the x-Lyapunov exponent for the system in (I^S)) is the time average of K{t) which is 
equal to —r(l+p<mo+p>TOi), where p< (p>) is the fraction of time that |a;s(t)| < 1 (|a;s(t)| > 1), where p< +p> = 1. 
Hence we have that the x- Lyapunov exponent is negative if p< < (1 + mi) /{mi — mp) = 2/3. From numerical solution 
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for the synchronized motion we find that this condition is indeed satisfied. On the other hand, the system (|26p is 
equivalent to s = — s — qs, where s — 5y(2), which converges toward the origin (0,0) with Lyapunov exponents, both 
equal to —1/2. Thus the synchronization of the network in (|23|) and p4)) . is ensured for sufficiently large networks. 




FIG. 6: The master stability function associated with the system (|27p as function of the real parameter A, varying between 
and 1. 



We now investigate synchronization stability for the systems ((23|) and ([24]) for a case of a real spectrum for Q, but 
without assuming large N . By linearizing the system in (|23p and (|24p about the synchronous evolution, we obtain for 
Eqs. ini): 



d 
di 



^ Sx ^ 



( 



K{t) \r 
A -1 -1 
-6.02 



&y{i) 



(27) 



In Fig. [21 we evaluate the master stability function associated with the system in (P7)) as function of A. The figure 
shows that, if all Ai in A' lie in the range (0, 0.7), synchronization will be stable. Moreover, since the master stability 
function (in Fig. [S]) becomes positive as A increases, the stability of the synchronous evolution depends only on 
^max] i-e., if M{Xmax) IS negative, the synchronous evolution is stable. Furthermore, we see that the large N limit is 
reasonably well satisfied for Xmax < 0.2 (i.e., M(A) at A = and at A < 0.2 are approximately the same). 

As a first example, we now consider a specific network where N is small. In particular we consider the network 
shown in Fig. [1] for which we have shown that the spectrum of Q is given by Eq. (|17p . Thus Xmax = \/w, and 
the spectrum of Q is real. Fig. [7] shows the synchronization error at large time as function of Xmax for w varying 
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between and 1 in steps of 0.01. We see that, in accord with our stabihty result from Fig. [6l stable synchronization 
of the chaotic motion is obtained if Xmax < 0.7. In obtaining Fig. [7l we initialize the variables Xi, y(i)i^ y{2)i randomly 
in a;i > on the synchronized chaotic attractor ([5]). For these initial conditions, we find that the time asymptotic 
synchronous motion is on the chaotic attractor of the system (rather than on one of the two fixed point attractors). 




FIG. 7: The network shown in Fig. [T]is modified as function of w. The plot shows the sum of the values of the errors Ex + Ey 
as function of the corresponding Xmax- 
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FIG. 8: The plot shows the transient evolutions Xi(t),i — l,..,Nx and yj(i,2)(i),J ~ i,---,Ny. The network parameters are as 
follows: Nx = 200 and iV„ = 300, p = 0.05. 



As a second example, Fig. |8] shows numerical results for a random network (case (i) of Sec. II. C) with Nx = 
200, Ny — 300, p — 0.05 corresponding to Xmax — 0.58, using the same type of initialization as in the previous 
example. Since Xmax < 0.7, Fig. [H] predicts stability, as is in fact seen in Fig. [51 In this figure the evolutions of all 
the randomly initialized systems is plotted versus time. It is seen that good synchronization is achieved by t > 10. 
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III. NETWORKS OF DISCRETE TIME SYSTEMS 

In this section we present a general analysis of two-group, bipartite network synchronization for discrete time 
systems. We assume the evolution of our discrete time network to be described by the following set of equations: 



(28) 



where Xi {yj) is a (riy) dimensional vector. Requiring Aij and Bij to satisfy conditions (|2ap and (j2bp . we see that 
multi-synchronous motion is possible and is described by the equations, 



yr' = G{y':) + L{x:). 

Linearization of (|28p about the synchronization manifold leads to: 



(29) 



(5x^+1 = DF{xs)5x^ + DH{y,) ^ A^.Sy"^ , i = 1, iV, 
Jy^'+i = DG{ys)5y'; + i?i(x,) £ B,,5x^, j = 1, Ny. 



(30) 



Similar to our previous analysis, we set = CxiSxn and (5j/" = Cy-Syn, where and Cy^- are appropriate scalar 
coefficients. Substitution of these into ([50)1 yields 

5x^+^ =DF{xs)5x'' + DH{y,)Y^2hL5yr^^ z = l,...,iV, (31a) 

5r+'=i?G(y,)Jr+^i(^.)E^^'^^"' J = l,-,^^y (31b) 

Then, following Sec. II, in order for piap (respectively (jSlbp ) to be satisfied for all i (respectively j) we require 
that (cxi)~^'Y^j AijCy. ~ {^yj)^^^i-^ji'^xi ~ where A is independent of both i and j. After defining the vector 
c — {cxi, Cx2, c-XN^ 1 Cyi, Cy^, Cy^^ ), the above conditions may be rewritten as Qc — Ac (as in Eq. 
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This lets us formulate the following master stability function problem 

= DF{xs)Sx'' + XDH{ys)Sy" 
Sr+' = DG{ys)Sy" + XDL{xs)Sx\ 



(32) 



As in Sec. II, the master stability function M(A) associates to each (possibly complex) the maximum Lyapunov 
exponent of the system in (|32p . The synchronous solution is stable if Af(A,;) < 0, for A^ in A'. 

IV. MORE GENERAL NETWORK TOPOLOGIES 



In this section we consider the case of more general network topologies. Specifically, we remove the constraint that 
the network is bipartite and we allow connections within the groups. We find that stable multi-synchronous evolutions 
are still possible and can be enhanced when intra-group connections are allowed. 

As an Example, we start by considering the following bipartite system. 



Ny 

ii{2) = 0.2 + a;,(2) (a;,(i) - 8.5), 



i = l,...,N^; 



(33) 



Vj = 0-2yj +^Sjia;,(i), j = l,...,Ny 

i=l 



(34) 



where A and B satisfy Eq s. Pa|) and (I2bp . In the synchronization manifold Eqs. ([55)) and ([M)) yield the following 
chaotic Rossler system {l8|. 



Xsil) = -Xs{2) - Vs, 

is{2) = 0.2 + a;s(2)(a;s(i) - 8.5), 
ys = 0.2ys + a;s(i). 



(35) 



Assume that the spectrum in A' includes zero as an eigenvalue. For such a case we see from (j34p that the y-component 
of the master stability equation pip yields dy — 0.2Sy, giving a Lyapunov exponent of 0.2 > 0. Thus the synchronized 
state is unstable for any network whose spectrum contains zero or small eigenvalues. 
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We now ask how this situation is affected by the presence of connections within a group. In order to illustrate this, 
we consider a case in which the system p3p . p4|) is modified by adding connections within the group Sy, leading to 
the following network equations, 

j=i (36) 
±02) = 0.2 + Xi(2)ixi(i) -8.5) i = 1,...,N^. 



i/j =0.2yj +J2BjiXi(i) - ayy^^CjkVk j = l,...,Ny, (37) 

1=1 fe=i 

where C = {Cjk} is a Laplacian matrix: "^^Cjk — for all j. It is important to note that the diffusive coupling 
term, CjkHk, is null in the synchronization manifold, where the dynamics is governed by the Rossler equations 
([55)1 . In what follows we consider a network with = 200, Ny = 300. We generate A and B randomly as described 
in Sec. II. C, case (i), with p^jy — Pyx = 0.10. We generate the Laplacian matrices randomly taking Cjk for j ^ k to 
be one with probability pyy and otherwise. Fig. [5] shows the effects of varying Uyy on the two quantities and Ey^ 
defined in Sec. II. We see that the network becomes synchronized for values of ayy > 0.1, indicating that diffusive 
intra-group coupling can be effective in enhancing the network synchronization. Note that, as Fig. [5] shows, though 
the diffusive terms are added only to the y-systems, synchronization applies for both the systems in Sx and Sy (in 
particular, what is observed is that both the systems synchronize in the chaotic Rossler evolution). However, we wish 
to emphasize that the master stability function approach presented in Sec. II, is inadequate for assessing the stability 
of the synchronous evolution when both intra-group and extra-group connections are allowed in the network. 

As another Example, we consider the system given by Eqs. (|23p and (|24p . introduced in Sec. II. The network 
topology is represented in Fig. [1] where here we consider the particular case oiw = 0.9. Thus, since Xmax — \/w > 0.7, 
according to the master stability function shown in Fig. [51 the network is not expected to synchronize. This is indeed 
what is shown in the left panels of Fig. (TUl where the systems trajectories of the x-nodes and y-nodes are shown to 
follow different evolutions. Now we consider whether it is possible to synchronize the network of these systems by 
adding diffusive couplings between systems in the same groups. Namely, we add a single bidirectional diffusive link 
between the two x-nodes and we assume a coupling constant equal to 2. That is, we add a term 2{x2 — xi) to the right 
hand side of Eq. for xi and a term 2{xi — X2) to the right hand side of Eq. for X2- As shown in the right 
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FIG. 9: The synchronization errors Ex and Ey versus (Jyy for Eqs. (j36|l and (|37|1 for a random network with Nx — 200 and 
Ny = 300, pxy = Pyx =0.1 and pyy = 0.15. 

panels of Fig. IIOI the network is now observed to synchronize on a multi-synchronous chaotic evolution. In particular, 
the equations for this evolution arc is = -r{xs + h{xs)) + ry^^i), 2/^(1) = -ys(i)+ys{2)+Xs, ys{2) = ~(iys{i) (again, 
observe that the diffusive coupling term is zero in the synchronization manifold). 



V. CONCLUSIONS 



Motivated by the common occurrence in applications of multi-synchronous motions in ensembles of interacting 



systems characterized by different dynamical behaviors [1|, [3, y, IJ, |5|j LZIi l9 



IC 



11 
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221, we have addressed 



the issue of how these systems can stabilize in distinct (possibly chaotic) synchronous evolutions. This form of 
synchronization is distinct from both diffusive coupling synchronization [16] and replacement synchronization 23|. 

By considering the underlying network of connections among the systems, we report conditions for the existence 
of a synchronization manifold. In the case of bipartite network topologies (i.e., when there are two communities and 
network links only connect nodes in different communities) we studied the stability of the synchronization manifold 
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FIG. 10: The left plots show the time evolutions Xi(t),i — 1, .., 2, and J/j(i,2) (i), J = 1, 3, for the bipartite network represented 
in Fig. [T] . The equations are those in (|23|l . (|24p . The state of the systems at the final time t* = 100 is shown by asterisks. 
The right plots show the state evolution of the network in the case when a bidirectional diffusive link with associated couphng 
constant equal to 2 is added between the two x-nodes. It is seen that the presence of the added link causes the network to 
synchronize. 



by means of a master stability function approach. In so doing, it was possible to decouple the effects of the network 
topology from those of the dynamics at the network nodes. We also presented an extension of our approach to discrete 
time systems. 

Finally, we considered examples of the case of more general network topologies, where links are also allowed to fall 
within each community, and we reported numerical evidence that the presence of diffusive couplings among nodes 
within the same community can enhance the network synchronizability 25!. 

We believe this paper represents only a first step in the study of multiple synchronization of complex networks. We 
hope that our work will stimulate further research efforts to address this issue in the future. 

This work was supported by ONR (Physics), NSF (PHY0456249), and by a MURI contract administered by ONR. 
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